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We confirm recent numerical results of echoing and mass scaling in the gravitational collapse 
of a spherical Yang-Mills field by constructing the critical solution and its perturbations as an 
eigenvalue problem. Because the field equations are not scale- invariant, the Yang-Mills critical 
solution is asymptotically, rather than exactly, self-similar, but the methods for dealing with discrete 
self-similarity developed for the real scalar field can be generalized. We find an echoing period 
A = 0.73784 ± 0.00002 and critical exponent for the black hole mass 7 = 0.1964 ± 0.0007. 
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J> 1 I. INTRODUCTION 

o : 

I Recently, Choptuik, Chmaj and Bizori [0 have studied the gravitational collapse of an SU(2) Yang-Mills field 
restricted to spherical symmetry. Near the boundary in initial data space between data which form black holes and 
data which do not ("critical collapse") they found two regions with qualitatively different behavior. In "region I" they 
found a mass gap, with the minimum black hole mass equal to the mass of the Bartnik-McKinnon solution, while in 
' "region II" they found the mass scaling and echoing which are by now familiar from critical collapse in other matter 
models. The two kinds of behavior are reminiscent of first and second order phase transitions. 

We review type I and type II behavior in section II. We shall see that each type of behavior can be understood 
through the presence of an intermediate attractor. The type I attractor is the well-known Bartnik-McKinnon solution, 
which is static and asymptotically flat. The type II attractor is self-similar and was not known before. After having 
\q ' derived the necessary field equations in section III, we construct it as an eigenvalue problem in section IV of this 
^\ . paper. In section V we construct its linear perturbations in another eigenvalue problem and verify that only one 
of them is growing. This allows us to calculate the critical exponent governing the mass scaling semi-analytically, 
q-( without numerical collapse simulations. In section VI we summarize our results, which are in good agreement with 
collapse simulations, discuss how the Einstein- Yang-Mills system differs from other systems in which critical collapse 
j^jrjl has previously been studied, and put the present paper into perspective. 

The analytic and numerical methods of this paper are a generalization of those developed for the spherical scalar 
field in In contrast to the scalar field system the field equations contain a length scale e _1 (in units c = G = 1), 

where e is the coupling constant in the Yang-Mills-covariant derivative D a = V a + eA a . The presence of a scale in 
the field equations excludes the existence of an exactly self-similar solution. Instead we make a series ansatz for a 
solution which becomes self-similar asymptotically on spacetime scales much smaller than e — 1 , or equivalently for 
curvatures much greater than e 2 . The echoing period A is determined by the leading term of the expansion alone. 
For the linearized equations we also make a series ansatz, but the spectrum {A^} of Lyapunov exponents is once more 
determined by the first term of that series alone. Moreover, to calculate the first term of the perturbation expansion 
one only needs to know the first term of the background expansion. Therefore the higher terms of either expansion 
are not required in order to calculate both the echoing period A and critical exponent 7 exactly, and will not be 
calculated here. 
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II. TYPE I AND TYPE II CRITICAL PHENOMENA 



Here we summarize the findings of Choptuik, Chmaj and Bizori |l| and explain them in dynamical systems terms. 
The purpose of this section is to show that critical phenomena are, in hindsight, easy to explain, to stress the 
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mathematical similarities between type I and type II critical behavior, and to motivate the more technical calculations 
in the following sections. 

For introductory reviews of critical collapse, see Very briefly, one wants to study the limit in phase space 

between initial data which eventually form a black hole and data which do not. Choptuik Q| pioneered the method 
of (numerically) evolving initial data taken from one-parameter families of initial data which cross this boundaries, 
families of data, in other words, which form a black hole for large values of the parameter, p, but not for small values. 
Generic families have this property. By bisection one can numerically determine the critical value p* of p for a given 
family. 

For the spherically symmetric massless scalar field Choptuik found that the black hole mass could be made arbitrarily 
small, and scaled like M ~ (p — p*) 7 , with 7 ~ 0.37 the same for all families of data. Furthermore, before forming the 
black hole, the time evolution from all data with \p — p*| sufficiently small, from all families, approached one universal 
solution with the strange property of being periodic in the logarithm of both r and t, or <f>(r, t) = <j)(e A r, e A t), with 
a period of A ~ 3.44. The smaller \p — p*|, the more "echos" were visible before the black hole formed or before the 
fields dispersed to infinity. 

For the spherically symmetric SU{2) Yang-Mills field, Choptuik, Chmaj and Bizon jj] found the same behavior, 
with A ~ 0.74 and 7 ~ 0.20 in some region of initial data space, and called it "type II behavior" because the black 
hole mass resembles the order parameter in a second-order phase transition. In another region of initial data space 
they found that as p went through p*, the black hole mass jumped to a finite value instead of showing power law 
behavior. Before the black hole formed, the solution approached a regular static solution of that mass, and remained 
there for the longer the smaller \p — p» was. They named this "type I behavior" , in analogy with a first-order phase 
transition. 

Without keeping up the sense of mystery any longer, we now explain these phenomena in dynamical systems terms 
through the presence of an "intermediate attractor" , a solution which has precisely one growing perturbation. 
We shall introduce a compact notation which focuses on the essential similarities and dissimilarities of the two types, 
while hiding the details, and which will also be useful later on: by Z we denote the vector of variables of some 
first-order form of the field equations, such that, for example, the complete field equations in spherical symmetry can 
be compactly written as F(Z, Z >r , Z, t ) = 0. The presentation is best begun with type I behavior. 

Type I behavior is dominated by the Bartnik-McKinnon solution [hoj , which is static, spherically symmetric and 
asymptotically flat. It has exactly one unstable perturbation mode which makes it an intermediate attractor in 

dynamical systems terms. Let Z^\r) denote this solution. As it depends only on r, its general linear perturbation 
must be of the form 5Z = Yli^Lo ^* eAi * $iZ(r), where the Ci are free constants. There is exactly one unstable mode, 

that is 9?Ai > and KAi < for i = 2, 3, Furthermore, it is known that the final state arising from initial data 

Zo(r, e) = Z*(r) + e 5\Z(r) is a black hole for one sign of e, and flat space with outgoing waves for the other. Let p 
be the parameter of a one-parameter-family of initial data, such that for p > p* a black hole forms, and for p < p» 
the solution disperses. Then for (p — p„) sufficiently small, the time evolution of data from the family enters an 
intermediate asymptotic regime of the form 

Z(r,t)^zi 1 \r) + ^-^) (p-p*) e^ h 8?Z(r), (1) 

where the decaying perturbations (i > 2) can already be neglected, and where we have approximated C\ — C\(jp) 
linearly. This solution leaves the intermediate asymptotic regime to form a black hole or disperse at a time T when the 
amplitude of the perturbations, dCi/dp (p — p*) e AlT , has reached some small fiducial value e. This gives a lifetime 

T ~ — Aln(p —p*)+c (2) 

of the meta-stable state, where c depends on the one-parameter family, but where A = l/A^ is universal. This was in 

fact observed by Choptuik, Chmaj and Bizon, and used to estimate Aj 1 , in good agreement with perturbation theory. 

The intermediate attractor in regime II is self-similar instead of static. It is in fact discretely self-similar, but that 
is not an essential detail, and for clarity of comparison with regime I we pretend in the introduction, and only here, 
that it is continuously self-similar. We also disregard the fact that the self-similarity is only asymptotic on small 
scales. In suitable coordinates Zi then depends only on r/t instead of only on r. Its general linear perturbation 
must be of the form SZ = J^So C 1 ; t Xi 8iZ(r). Once more there is exactly one growing mode. (This is required for 
the explanation we are about to give to work, and we will demonstrate it explicitly in section IV.) The intermediate 
asymptotic regime for type II behavior is 

Z(r,t) ~ Zi u \r/t) + ^-(p*) (p-p.) t A ™ S^Z(r/t). (3) 
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One now argues from scale-invariance |3|]9|,|l2|,|3|] that the black-hole mass M is proportional to T, and obtains for the 
black hole mass 

InM =7ln(p-p») +c, (4) 

where 7 = — 1/A^ is universal, and c is a family-dependent constant. 

If the scale-invariance is only asymptotic, as it is for scalar electrodynamics or Einstein- Yang-Mills, the scaling 
argument to calculate the black hole mass goes through unchanged If the critical solution is discretely self- 

similar, as for the scalar field or the model considered here, with an echoing period of A in the logarithm of the 
length and time scales, a periodic "wiggle" Q or "fine structure" |1J] is superimposed on the mass scaling law, which 
becomes 

InM = 7ln(p-p,) + c + *[ln(p-p«) + c/7], (5) 

where $ is a universal periodic function with period A/ (27). (Note that the same family-dependent constant c appears 
a second time in the argument of The form of the critical solution and its perturbations is also more complicated, 
and will be discussed in section IV.A and V.A respectively. 



III. FIELD EQUATIONS AND SCALING VARIABLES 



In this section we introduce coordinates and field variables for the spherically symmetric Einstein- Yang- Mills system 
that are adapted to type II behavior, where scale-invariance plays a crucial role. In the following we consider only 
type II behavior, and no longer write the index (II). We adopt the conventions and notation of M, which include 
making both the Yang-Mills field and the coordinates r and t dimensionless by absorbing suitable factors of G, c and 
e into them. 

The spherically symmetric spacetime metric is written as 

ds 2 = -a 2 dt 2 + a 2 dr 2 + r 2 (d$ 2 + sin 2 d dip 2 ) , (6) 

where a and a depend only on r and t, and the Yang-Mills field strength is given in terms of a single potential W (r, t) 

as 

F = dW A (n dd + r 2 sintfd^) - (1 - W 2 )r 3 d^ A sintidy, (7) 
where n are the Pauli matrices. In order to write the field equations in first order form, we define 

$ = VF r , 11= -W t . (8) 
a 

The complete field equations, reduced to spherical symmetry, are 

r$ )t = rf-n) , (9) 

V a / ,r 

rU t =r(-^) +aar- 1 W(l - W 2 ), (10) 

V a / ,r 

r^- = 1(1 - a 2 ) + $ 2 + IT 2 + ia 2 r- 2 (l - W 2 ) 2 , (11) 
a 2 2 

r ^L = I( a 2 _ 1) + $2 + n 2 _ I a 2 r -2 (1 _ ^2)2 (12) 

a 2 2 

r^l = 211$. (13) 
a 

These equations are the Yang-Mills equation, and three of the four algebraically independent components of the 
Einstein equations. The fourth component is obtained by combining derivatives of the other three, and is therefore 
redundant. 

In order to construct a discretely self-similar solutions, we follow in defining new coordinates 

r = ln(-t), C = ln(-^) -&(r), (14) 
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where £o is a periodic function to be determined, with period A. (This definition differs slightly from || in that t and 
r are dimensionless, and that t is negative.) The resulting spacetime metric is 



ds 1 



T |-a 2 dr 2 + e 2(c+ « o) a 2 (d( + (1 + Q drf + M 2 + sin 2 <& dip 2 } 



where a and a are now functions of £ and r, and where £o = £o{t) and £q = d^o/dr. As discussed in ||, 
self-similarity is equivalent to a and a being periodic in r. In the field equations we make the replacements 



d_ 

dr 



d_ 



r H. - P C+Co(r) 
dt ~ 



r = e ^+C+Jo(r) 



(15) 
discrete 

(16) 



to transform to the new coordinates. 

We shall be looking for a solution in which a and a are periodic. What does this mean for the matter variables <E>, 
II and Wl The Einstein equations suggest that $ and II should be periodic too, but W cannot be periodic because 
of the explicit presence of the factors e T in the equations. Nor can we simply absorb such a factor into the definition 
of W to make it periodic. This means that the equations have no nontrivial self-similar (periodic) solutions. The 
physical reason is the presence of the length scale e _1 in the problem, which is only hidden by the dimensionless 
variables. Following a suggestion by Matt Choptuik Ea], we define a new scalar field S by 



W = l-rS. 

With this definition, the two potential terms arising in the field equations, 

r -1 W(l - W 2 ) = (1 - rS)(2 - rS)S, r~ 2 (l - W 2 ) 2 = (2 - rS) 2 S 2 , 



(17) 



(18) 



split into the sum of a term which no longer contains r explicitly, plus terms containing positive powers of r, which 
become negligible on small spacetime scales (as r — > or t — > -co). Two further definitions, namely II± = II ± <J> 
and g = a/a, will be useful because g alone determines the ingoing and outgoing null geodesies, and II + and II_ are 
the components of the matter field propagating along them. 

In the following, we use the coordinates C and r, and the fields Z = {a, g, II+, II_, S}. In these variables, the 
complete field equations, including the definitions of II + and II_ in terms of S, are 



n±,c 
s,c 





Te^+« 5 n ±:r + Cn± T a 2 (l - e T +«+«°5)(2 - e T +^°S)S 



lT(l + S)e f+6 9 



a 
2 
Cg 



{c + ii 2 



n 2 



v 2 ; " 



-(C+Co) 



n_), 

g -^(Ul-Ut)-(l + e )Uc + Ul+U 2 



Q = S, 



-(C+€o) 



- 1 -(n + + n_) + (i + ^) 



^ + -(n + -n_) 



where 



C = l-a 2 +a 2 (2-e 



(19) 

(20) 
(21) 
(22) 
(23) 

(24) 
(25) 



As suggested by the way we have written the equations, the first five can be treated as evolution equations in £, 
with periodic boundary conditions in r, and the last two as constraints which are propagated by the evolution 
equations. Note that now only positive powers of e T appear explicitly, so that in the limit r — > -co we are left with 
a set of nontrivial, scale-invariant equations for Z. The terms multiplied by e r are "irrelevant" in the language of 
renormalisation group theory |l6| ] . 

The equations are invariant under W — > — W, and the potential for W has the two minima W = ±1. In ([IT]) we 
have assumed that W — » 1 asymptotically. A solution tending to W = — 1 can be trivially obtained from one tending 
to W = 1 by changing the sign of W, $ and II, while leaving S, a and a unchanged. The field equations are left 
unchanged. 
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IV. BACKGROUND SOLUTION 



A. The eigenvalue problem 



In this section we construct the solution zi U \r, t) which dominates type II behavior. To find a solution which is 



asymptotically self-similar in the limit r 
field equations, we make the ansatz 



-co, that is on spacetime scales smaller than the intrinsic scale of the 



n=0 



e BT Z, n (C,r), 



(26) 



where each Z m is periodic in r. Z*q is the solution of a nonlinear eigenvalue problem, with eigenvalue A, and boundary 
conditions arising from certain regularity requirements. Z*i is the solution of an inhomogeneous nonlinear boundary 
value problem, with source terms depending on Z»o- Similar boundary value problems completely determine all higher 
Z m recursively. 

In the following we are interested only in the equations for Z*q, and from now on we suppress the suffix *o on the 
components of Z*o, denoting n + »o simply by II + etc. (In the compact formal notation Z*q we keep the suffix.) The 
equations for Z„q are derived from those for Z above by setting the factor e T equal to zero at each explicit occurrence. 
We choose to evolve only 11+ , II _ and g in £ with eqns. (|l9| ) and (21), and to determine a and 5 at each new value 
of £ from the constraints, eqns. (E3|) and (^4|). The final set of equations, those we have solved numerically, is, in the 
simplified notation, 



n±, c 



9X = 



Te c+ ^. 9 n ±!T + (1 - of + 4a 2 S 2 ) U ± T 2a 2 S 

1 T (1 + &)e<+«° 5 
(l - a 2 + 4a 2 S 2 ) g, 



= a T + e-^g- 1 ^(U 2 



U 2 



(l+Q^(l-a 2 +4a 2 S 2 + n 2 + +Il 2 _) 



o = s, T - e-x+^g- 1 ^ (H+ + n_) + (1 + & 



s + -(n + -n_ 



(27) 

(28) 
(29) 

(30) 



All fields areperiodic in r with a period A that is to be determined as an eigenvalue. Here as in the example of the 
scalar field Q, the field equations are complemented by regularity conditions at the center r = (for t < 0), and 
at the past self-similarity horizon (the past light cone of the point (r = 0, t = 0), or r ~ —t). One can solve these 
boundary conditions in terms of free parameters. 

To make r = <^> ( = -co a regular center, we impose a — 1 and g = 1 there. We expand in powers of e^, and 
notice that a, g and II are even in that expansion (because they are even in r at r — 0), while S and $ are odd. We 
label the orders of this expansion by a suffix in round brackets to distinguish them from the orders in the expansion 
(p6|). The expansion coefficients can be given recursively in terms of one free periodic function 5(1) (r). To order e 3 ^ 
they are, giving II and $ instead of n + and II_ , 



a (0)(T 

.9(0) (t 
n ( o)(r 

0(2) (t 

.9(2) (r 

n (2) (r 



(3)1 



1, 
1, 

o, 

free. 

-25(i) 
25( 2 i ); 
0, 



5(1) - (1 + £ )5(i) 



f 

10 

-45( 3 ). 



{ e «° [n' (2) -2(i + e )n (2) 



85 ( 3 d} , 



These expressions are used to impose the asymptotic boundary condition at C 

C = Cleft- 



(31) 
(32) 
(33) 
(34) 
(35) 

(36) 
(37) 

(38) 

(39) 
(40) 

-oo at some small value of £ say 
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We use the remaining coordinate freedom, the choice of £o ( T ) , to move the self-similarity horizon to the coordinate 
surface £ = through the coordinate condition 

[l-(l+Qe^g] c=Q = 0, (41) 

which means that £ = is null, and impose analyticity there by the condition 

[-e^°gIl +tT + (1 - a 2 + 4a 2 S 2 )U+ - 2a 2 S] Q=Q = 0. (42) 

(This is a regular and sufficient condition, by the same argument already used in B.) 

These two constraints can be solved recursively after expanding, this time in powers of £. We denote the components 
of this expansion also by subscripts in round brackets. The two free parameters here are the periodic functions <7(o)(t) 
and IT_( )(t). From (f4l|), one obtains the algebraic identity 

g {0) = [e^(l + Q}-\ (43) 

To obtain the leading order coefficients of the other fields, we substitute ( ff3| ) into eqns. (|30|), (p9|), and (^7|) (upper 
sign), obtain 



o'(o) - (1 + Co) (n^( ) + I) fl(o) - (1 + Q (25, »;;,, = «>. ur.) 



S(o) + (1 + Q S (o) - (1 + ^o)n-(o) = 0, (44) 

'in, - 2 J "(0, 

n ' + (o) - (1 + £0) (l - 4) + 4a (o) 5 (o)) n +(o) + (1 + QHo) S (o) = 0» ( 46 ) 

and consider these as linear ODEs for SVqJj (Q(o)) _2 i an d H+(o) respectively. 

As in the scalar field case, we make the assumption that the metric variables a and g contain only even frequencies 
in r, and the matter variables II_ and S only odd frequencies. This is compatible with the equations for Z*q, 
but not with the equations for the general Z. If this symmetry did not hold, the right-hand side of eqn. would 
contain even terms in r, and among them generically a term constant in t. Then S would not be periodic in t, but 
would have a term linear in r, and through the Einstein equations this would be in contradiction to the periodicity 
of a and g, and hence the self-similarity of Z»o- 

The equivalent of the field S here is the scalar field 4> in the scalar field model, and for a massive or self-interacting 
a similar argument holds. The equations for a massless <f>, however, do not contain <f> itself but only its derivatives 
IT±. Therefore a linear dependence of cj) on r would not clash with spacetime self-similarity. Such solutions exist, and 
have been investigated in |[7| , but surprisingly the critical solution for the massless field is not of this kind, and the 
massless and massive (or self-interacting) scalar field are therefore in the same universality class. 



B. Numerical construction 

Our numerical method has been described in detail elsewhere By decomposing all fields in Fourier components 
with respect to r, the PDEs in r and ( go over into a (large) system of ODEs in the variable C for the Fourier 
components. ODEs in r alone, in the boundary conditions and the constraints, go over into algebraic equations which 
can be solved in closed form. A now appears as a parameter in the Fourier transformation of the r-derivatives. 

A solution of the field equations and boundary conditions exists only for isolated values of A, and we have found 
precisely one. The convergence radius of our relaxation algorithm is smaller than for the scalar field, probably because 
of the shorter period A, and instead of an ad- hoc initial guess we had to use collapse data kindly provided by Matt 
Choptuik []l5| to obtain a good enough starting value for the relaxation algorithm. 

We find good agreement of Z„q with the Z of a critical collapse simulation for —3.00 < r < —2.22 ]Fq| , which is not 
very surprising as we started our numerical search with these data, but nevertheless confirms that the ansatz (|26| ) for 
Z* is consistent and converges for small enough r, with Z*q the dominant term. 

To obtain error bars on the solution, we have checked convergence with the numerical parameters, Cleft; the number 
N of Fourier components and the grid spacing A£, by varying one of them at a time. 

Fig. [l] demonstrates quartic convergence with exp£i e ft, as expected from our expansion to order 0(exp 3£i e ft). This 
convergence breaks down at very small values of exp Cleft j due to the fact that all fields become very small. 

Fig. ^ demonstrates quadratic convergence with grid spacing in as expected from centered differencing of the 
(■-derivatives. This convergence breaks down at very small values of A£, probably because grid points get very close 
to the regular singular point £ = 0. 
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Convergence with N is rapid: The difference between results for N — 64 and N = 128 is already of order 10~ 6 . 
N = 64 is surprisingly small, given that it means only 16 odd Fourier components each to represent 11+ (r) and II_(t) 
and 16 even components for g(r) and 15 for £o( r )- (The component cos(47r/A) of £o( r ) is taken to be zero to fix the 
translation invariance in r of the equations for Z*q.) 

For the production run we have chosen ^ieft = —6.4, A£ = (1/80) (that is, 513 grid points) and N = 128. The 
solution Z*o(Ci r ) nas an estimated maximal error of ±2.3 • 10~ 4 and root-mean-square error of ±3.6 • 10~ 5 , in the 
region —6.4 < £ < 0. We obtain A = 0.73784 ± 0.00002. All three error estimates are dominated by the error from 
finite differencing in £, with the estimated error from expanding around £ = — oo somewhat smaller, and the error 
from using a finite number of Fourier components in r much smaller. 



V. LINEAR PERTURBATIONS AND CRITICAL EXPONENT 



A. The eigenvalue problem 

In this section we construct the one linear perturbation of the critical solution that grows with decreasing spacetime 
scale, as r — * — oo, with the purpose of calculating the critical exponent for the black hole mass in critical collapse. 
The evolution equations for a linear perturbation SZ of the background critical solution Z* are of the general form 

SZ X = A 5Z, T + (B + e T C + e 2r D) SZ. (47) 

The constraints are of the same general form, but with the left-hand side equal to zero, and the following considerations 
apply equally to them as well. The perturbation equations differ from those for the scalar field model through the 
explicit appearance of e T in the equations, and the fact that the coefficients A, B, C and D are not periodic in r 
because the background solution Z* is not. Like Z#, the coefficients A, B, C and D derived from it admit an expansion 
of the form 

oo 

A = J2e nT A n (C,r), (48) 
n=0 

where the A n are periodic. In this expansion, the leading terms Aq, Bq, etc. depend only on the leading term Z*q of 
the background expansion. 

As for the scalar field model, we make the ansatz || 

oo 

6Z(t,T) = J2 C i e XiT 6iZ(C,T), (49) 

2 = 1 

where the Ci are free coefficients, and the A,; are a discrete set of complex numbers, which arc determined as eigenvalues 
of a new, linear boundary value problem. Clearly the 5iZ obey the equations 

5iZ x = A 8iZ tT + (B + X t A + e T C + e 2r D) 6 t Z. (50) 

In the massless scalar field model, the 5iZ could be assumed to be periodic in r. In the presence of a scale, this is no 
longer possible, and we have to expand each SiZ once more as [ p"3| 

oo 

S i Z(C,T) = Y f e nT S in Z(C,r), (51) 

71=0 

where only the individual coefficients 5i n Z are periodic. This expansion is exactly analogous to (p6|). The Si n Z obey 
a coupled set of equations which can be derived from ( |50| ) in a straightforward bookkeeping exercise, after inserting 
the expansion (|48|). These equations are complemented by regularity conditions at £ = — oo and C — 0. The equations 
for the SioZ are simply 

SioZx = A Q 6iaZ tT + (B + XiA ) S i0 Z. (52) 

This equation, together with the boundary conditions, already determines the spectrum {A;}. The other Si n Z obey 
inhomogeneous equations and can be determined recursively, but here we are interested only in the spectrum. This 
also means that we only need the leading term Z* of the background expansion. 
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Writing down the field equations ( |52| ) for the SioZ is straightforward. As we have seen, one simply linearizes eqns. 
( ^7| - |30| ) for Z*o, and then makes the replacement 5^*o.r — ► 5ioZ, r + A; 5ioZ, T , which follows from the definition (|is|). 
Writing a for et*o, etc., and 5a for <5ioa, etc., to keep the notation simple, we obtain 



oTI ± , c = 



[1 T (1 + £o)e C+e °] 1 {Te C+Co [II ± , r <5 9 + g(6TL±, r + A, SIL±)] + (1 - a 2 + 4a 2 S 2 ) 5T1 ± 
-f 2a(45 2 - l)n± 5a + 8a 2 STI ± SS T 2a 2 5S T 4a5 6a ± (1 + £ )e C+?o n ±iC <^}> 



<5g, c = (1 - a 2 + 4a 2 5 2 ) <5.g + 2a(45 2 - l).g (5a + 8a 2 Sg <55, 
= 5a T + 

z 

1 



(53) 
(54) 



A, + i e -(C+«o) 5 - 1 (n 2 , - n 2 _) - 1(1 + Q(l - 3a 2 + 12a 2 5 2 + Ii\ + n 2 _) 



<5a 



+ ) e -(C+€o) 



-|(n 2 - n 2 _).g- 2 a <5g + ff~ 1 a(II + <5II + - II oTI_) 



-(1 + Co) [4a 3 S 5S + a(n+ (5n+ + II_ oTI_)] j , (55) 
= <55 r + (A. t + l+^) ^+l{e-«+«°) [ 5 - 2 (n++n_) (5.g-«r 1 (<5n + +<5n_)] + (1 + £ )(<5II + - <5II_)} . (56) 



Similarly, we obtain the expansion around £ = — oo of the S^Z by linearizing eqns. ( |3l| - |4C| ) and then making the same 
replacement, at each order in e^. The nonvanishing expansion coefficients to 0(e 3 ^) are: 



SS(i) (r 
(5$ ( i)(r 
(5a (2 )(T 

OT (2) (t 

<55( 3) (t 
<5$ (3) (t 



free, 

-2 (55(i), 
45(i) <55 ( i) 
P 4o 



«55( 1) + (A-l-e ) 
i {e«° [«5n' (2) + (A - 2 - 2£ ) <5II (2) <5] + 245 2 1} d5 (1) } , 
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-4 (55(3) 



(57) 
(58) 
(59) 

(60) 

(61) 
(62) 



As the linearized regularity condition at £ = we impose the vanishing of the numerator of (^3[ upper sign) . There 
is no linearized equivalent of the coordinate condition (^lj) , as we have fixed the coordinate system already when the 
background was calculated. (In other words, C = is null only on the background spacetime, not on the perturbed 
spacetime.) The one boundary condition at £ = can be solved recursively in terms of two free periodic functions 
$9(o)(t) and i5n_ (0) (r), from 



«55'+(A, + l + £ ) 55(0) + (1 + Q 



-<5II_ (0) + -g 1 (n +(0) + II_(o)) 5g {0) 



= 0, 



5a', 



(0) 



Ai + (1 + Q ( -n 2 _ (0) - - + -a 2 (Q) - 6a 2 0) 5 ( 2 0) 



(5 rt 



(0) 



n' 



7^(o) a (o) ( n +(o) - n -(o)) s 9(0) +4a^ 0) 5 (0) 5S (0) + 2a (0 )II_ { o) oTI_ (0 ) 



= 0, 



A, - (1 + £o) (l - a 2 0) + 4a 2 )5 2 ))] 6IL +{0) + g (0 ) [n' +(0) - (1 + £ )n +(1) ] 5g {0) 
+ (1 + Q {2a (0 ) [(l - 45 2 0) ) n +(0) + 25 (0) ] ^(o) + 2a 2 0) [l - 45 (0) n +(0) ] «55 (0) } = 0. 



(63) 



(64) 



(65) 



The suffix (0) denotes the leading term in an expression in powers of £ around ( = 0. We still need to calculate the 
background term (IL^f)/^ = II_|_(i) in eqn. (J65|). To do this, we expand eqns. @ lower sign), (§|) and (|| to O(C), 
evaluate the resulting algebraic expressions 



n 



■9(1) = C \o)9(o), where C (0) = 1 - a 2 0) + 4a 2 0) 5 ( 2 0) , 
-(i) = 9 C 1 +^o) _ln '-(o) + CoII-o + 2ft( 0) 5( ) , 



(66) 
(67) 
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S(i) - -5 (0) - - (n +(0) 



(0)J 



O(l) 



a (Q) ^ 



(0) 



n 



(0) 



n 



"(0) 



(68) 
(69) 



(alternatively, we could have obtained Sm and am from expanding the constraints (55) and (|56j)), and finally solve 
the linear ODE 



nV (1) -(i + £o)(i + 2C ( o))n +(1) 



(i + c (0) )nv (0) 



-(l + Q {2a (0) [(l - 4S ( 2 0) ) n +(0) + 25 (0) ] «(i) + 2af 0) t 1 " 45 (o) n +(o)] S (1) } = 



(70) 
(71) 



for n 



+ (!)• 



Linear perturbations which have the same r-symmetry as the background Z*q [S and II± odd frequencies, a and 
g even frequencies) decouple from those with the opposite symmetry. We call them even and odd perturbations 
respectively, and can treat them separately in the numerical calculation of the spectrum {A;}. 



B. Numerical construction 



Our numerical method is the same as in j3||. We evolve a basis of all linear perturbations compatible with the 
constraints at either one of the boundaries to a matching point, and look for zeros of the determinant of that basis as 
a function of A. A zero indicates the existence of a perturbation consistent with both sets of boundary conditions for 
that value of A. We have implemented this algorithm for both real and complex A. We have checked our results, for 
real A and even perturbations, with a relaxation algorithm, which is partially independent numerically, and in which 
A figures as an additional variable, which is balanced by fixing the perturbations as an additional boundary condition. 
The determinant in question is in fact a holomorphic function of A (because the field equations are real) , and this can 
be used to find its zeros and poles efficiently. 

We expect certain zeros and poles in the A-plane from the following considerations. Z*o is scale- invariant, and 
therefore invariant under the infinitesimal transformation 



Z* (r, t) -> Z* [(l + e)r, (1 + e)t] ~ Z* ((, r + e) ~ Z„ (C, r) + eZ^ T . (72) 

This corresponds to a gauge linear perturbation mode with Ai = and SiZ = Z*o,t- -Z*o is also invariant under time 
translation, 

Z*o(r, t) -> Z* (r, t + e) ~ Z* ((, r) + e e~ T [(1 + £ )-^*o,C ~ Z *0,t) , (73) 

corresponding to a gauge mode with Ai = — 1. Both gauge modes are even according to our classification. 

The ODEs, eqns. ( p3| , |64] , |65| ), are all of the form /' + gf + h = 0, where / stands for <5SVo), Sa^ and <$n + ( ) 
respectively. In all three equations g depends only on the background solution and is even, while h is linear in the 
perturbations, and has the same T-symmetry as /. It can be shown || that this type of equation has no solution 
when the average value (in r) of the coefficient g vanishes. As g in each case is of the form A + (background fields), 
this corresponds to a simple pole in the A-plane. These poles are not just due to the breakdown of a particular 
numerical method but indicate that for these values of A no perturbations exist which obey the boundary condition at 
C = 0. The poles arise only when the inhomogeneous term h, and in consequence the unknown /, have a nonvanishing 
average, that is when they are even. 

Calculating the average value of g for each of the three equations, we find that they vanish for A = — 1, A = — 1 — A 
and A = —A respectively, where A is the average value of 2(1 +£ )ri'5, s, with numerical value A ~ 0.1726. (We have 
used the fact that 




has vanishing average value, as it is the derivative of a periodic function, to simplify the averages.) In summary, for 
even perturbations we expect zeros at A = and A = — 1 (gauge modes), one more zero on the negative real line (the 
unstable mode) and a pole at A ~ —1.17. For odd perturbations we expect poles at A = — 1 and A ~ —0.17. 

The numerical calculation of the perturbation determinant as a function of A largely confirms the predictions: For 
even perturbations, on the negative real line we find a zero at Ai ~ —5.0, corresponding to the expected physical 
unstable mode, with a critical exponent of — 1/Ai ~ 0.2, as found in collapse simulations. We also find the expected 







zero at A2 = 0. We have verified that the corresponding 5iZ oc Z*o,t to high precision. We find the expected pole 
at A ~ —1.17, but accompanied by a zero very close by. For odd perturbations, on the negative real line we find the 
expected pole at A ~ —0.17. 

At A = — 1, for both even and odd perturbations, we do not find the expected zero and pole respectively because 
of a numerical problem which is discussed in the appendix. It does not affect our calculation of the perturbation 
determinant for values of A not close to — 1. The unstable mode at A ~ —5.0 and gauge mode at A = are clear 
enough, and we can use their convergence properties to obtain an estimate of the numerical error. 

Table 1 gives the values of A for the unstable mode Ai and the scale change gauge mode A2 as a function of the 
step size A£. The deviation of the numerical value of A2 from zero serves as one estimate of numerical error. It is 
larger than the other estimate, from the convergence of Ai, and we therefore adopt it as our definitive error estimate 
for Ax- We obtain Ai = —5.091 ± 0.017, from which we obtain for the critical exponent 7 = — 1/Aj = 0.1964 ± 0.0007. 



VI. CONCLUSIONS 



We have obtained the asymptotic form of the type II critical solution of Einstein- Yang-Mills collapse, its echoing 
period, and the critical exponent for the black hole mass, in a calculation similar to the one we made for the 
massless scalar field ^,|). The major new feature is the presence of the length scale in the Einstein- Yang-Mills field 
equations. In consequence, the critical solution and its linear perturbations are no longer self-similar, but become so 
only asymptotically on spacetime scales much smaller than the length scale of the field equations (on the order of 
the Bartnik-McKinnon EQ] mass). Here we have only calculated the leading term in the asymptotic expansions for 
the critical solution and its perturbations, but this is sufficient to calculate both the echoing period A and critical 
exponent 7 exactly. We find A = 0.73784 ± 0.00002 and 7 = 0.1964 ± 0.0007, while Choptuik, Chmaj and Bizon 
JjJ find A ~ 0.74 and 7 ~ 0.20 in collapse simulations. Data files of the background and unstable m ode from the 
production run are available through the WWW address http://www.aei-potsdam.mpg.dc/^gundlach. 

In the formalism we have developed here to deal with the presence of a length scale in the equations, the leading 
perturbation term, ^2 Cie XiT S^Z ', obeys field equations which are the linearized version of the equations for the leading 
background term, Z*o- Both sets of equations consistently describe a new physical system which is scale-invariant, 
and which is obtained from the original, scale-dependent, model in the limit where all fields vary on spacetime 
scales much smaller than the intrinsic scale of the field equations. In the language of renormalisation group theory, 
these equations are the short-scale fixed point of a renormalisation group transformation acting on the original field 
equations. In the case of a massless or self- interacting scalar field this fixed point is the massless scalar field |J. For 
scalar electrodynamics, the fixed point is the massless complex scalar without electromagnetism Jj"3"| ]. In both cases 
the field equations at the fixed point can naturally be associated with the Lagrangian of a simpler physical system, 
and the renormalisation group acts naturally on the dimensionful coupling constant. While the latter is still true 
for the spherically symmetric Einstein- Yang-Mills system, the equations (27-3CJ) are not the spherical reduction of 
some set of covariant field equations. The reduction to spherical symmetry does not commute with the action of the 
renormalisation group. 

We believe that in the present paper we have developed the most general formalism that will be be required to 
deal with critical collapse restricted to spherical symmetry, in allowing for self-similarity and the presence of a length 
scale. The generalization to more than one scale is trivial: the various scales can be written as a single scale times 
dimensionless numbers. The general formalism has already given rise to a bit of new physics: the calculation of critical 
exponents not only for the black hole mass but also its charge in critical collapse of scalar electrodynamics (Df ]. 

The example of Einstein- Yang-Mills collapse shows that one does not need scale-invariance of the field equations to 
have type II critical phenomena with the famous relation M ~ (j> — p*) 7 . Rather they can be found in some region of 
initial data space for any system. For astrophysical matter, these initial data are simply not realized in astrophysical 
collapse. 

Most remaining questions in critical phenomena go beyond the restriction to spherical symmetry. Do the spherical 
critical solutions found so far act also as critical solutions for generic, non-spherical, initial data? What is the angular 
momentum of the black hole formed from data with angular momentum in the limit where the black hole mass is 
fine-tuned to zero? Are there qualitatively new phenomena away from spherical symmetry? 



ACKNOWLEDGMENTS 



I would like to thank Piotr Bizon and Matt Choptuik for helpful discussions, and Matt Choptuik for making his 
data from collapse simulations available. This work was supported by a scholarship of the Ministry of Education and 



10 



Science (Spain). 



APPENDIX A: NUMERICAL PROBLEMS AT A ~ -1 



Calculating the determinant of even perturbations as a function of A, we do not find the expected zero but a pole at 
A = — 1, with an alternating series of poles and zeros accumulating towards —1 from below. There are no further poles 
or zeros immediately above. The positions of the poles are the same in the real and complex algorithms, but depend 
on the values of the numerical parameters N and 5£. A qualitatively similar picture arises for odd perturbations. 
These features can be explained as a numerical artifact as follows. 

We have checked explicitly that the gauge mode ([73]) obeys the constraint (p5|). When we try to reconstruct 5a of 
this mode from the constraint, however, the numerical result blows up at small £. To understand this, consider the 
equation 5a' + g 8a + h = 0, with g and h defined by eqn. (|55|). The Fourier algorithm that we use to solve this for 5a 
at each £ needs to divide the average of h by the average of g. As £ — > — oo, the average of g over r as a function of A 
and £ is A + 1 + 0(e 2< »), where the last term is positive. As A — > — 1 from below, this goes through zero at some small 
value of (. In the exact perturbation mode (|7^), the average value of h vanishes at the same rate with £ as that of g, 
but with small numerical errors this cancellation fails, and small numerical errors are magnified. In the calculation of 
the perturbation determinant this results in the observed, essentially random behavior for A J$ — 1. For A > — 1 the 
problem does not arise, as then the average value of g does not vanish for any £. 

We have not found a simple way of fixing this problem, as our algorithm relies in an essential way on reconstructing 
S and a and 5S and 5a from the constraints at each £. It does not affect numerical results however unless where the 
average values of both the coefficients g and h are very small, that is for A < — 1. (If only the average value of g is 
small, the resulting blowup in the perturbations is physical, as in the other poles we have discussed.) Calculating the 
perturbation determinant is not goal in itself, but only a means of finding the spectrum of linear perturbations. With 
the present method we can say with confidence that there is a zero at A ~ —5.0, and no other zeros for negative real 
A, apart from the two gauge modes. We could in principle be missing a zero (physical growing mode) at A < — 1, 
where the code is unreliable, and have to rely on evidence from collapse simulations that there is only one unstable 
mode. 
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FIG. 1. Quartic convergence of Z,o with exp£i e ft. Assuming that the error is E ~ Aexp4(\ c { t , the difference between two 
numerical solutions obtained with (j e f t and (j c f t + A£i e f t is E' ~ 4AA(i c f t exp4£i e f t . Therefore we plot here £'/(4A£i e f t ) ~ E, 
against exp£i e f t . Circles denote the maximal error, over all gridpoints and Fourier components, squares the root mean square 
error, and triangles the error in A. N = 64 and A£ = 0.1. The production value £i e ft = —6.4 corresponds to the fifth point 
from the left. 
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FIG. 2. Quadratic convergence of Z*o with A£. As a measure of the error at A(, we compare with A£/2. Circles denote the 
maximal error, over all gridpoints and Fourier components, squares the root mean square error, and triangles the error in A. 
N = 64 and £i e f t = —6.4. The production value A£ = (1/80) corresponds to the second point from the left. 



TABLE I. Convergence of A with step size in £. Ai is the Lyapunov exponent of the one growing mode. Its negative inverse 
is the critical exponent 7. A2 is the exponent of the scale change (r-translation) gauge mode. It must be zero and serves as a 
check on the numerical error. Note that the numerical grid (and number of steps) is the same for the background as for the 
perturbations in each case. The range of £ is —6.4 < £ < in each case. 

Number of steps Ai A2 

32 -5.1318584162589 0.13720816860828 

64 -5.0828194924873 -1.4932102080912E-02 

128 -5.0891816495598 -1.1519697001693E-02 

256 -5.0910562847918 4.7061190684163E-03 

512 -5.0913625725286 1.6758611037677E-02 
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